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Abstract 

The  objective  of  the  project  is  to  develop  and  improve  upon  three-dimensional  (3-D)  seismic  velocity  mod¬ 
els  of  the  Earth  and  to  utilize  such  models  for  improving  the  locations  of  events  recorded  at  regional  and 
teleseismic  distances.  Our  previous  experiments  performed  using  3-D  models  of  varying  resolution  show 
that  the  accuracy  of  event  location  with  respect  to  ground  truth  does  not  necessarily  improve  with  an  in¬ 
creasing  number  of  free  parameters  in  a  model.  One  reason  for  this  may  be  the  lack  of  waveform  data  used 
for  constructing  recent  high-resolution  models,  since  the  maximum  sensitivity  of  seismic  travel  times  to 
structure  is  in  the  middle  to  lower  mantle.  In  this  paper  we  develop  a  new  joint,  P  and  S  velocity  model 
parameterized  in  terms  of  radial  and  horizontal  cubic  splines,  using  a  combination  of  direct  and  differential 
travel  times  and  surface  wave  phase  measurements.  The  accuracy  of  the  travel  time  data  has  been  improved 
by  relocating  the  events  using  a  previously  existing  3-D  mantle  P  wave  model.  The  cubic  spline  param¬ 
eterization  conveniently  allows  more  detailed  models  to  be  inserted  within  it  for  use  of  regional  phases  in 
event  location.  Our  preliminary  model  shows  substantial  differences  in  compressional  velocity  in  the  upper¬ 
most  mantle  compared  with  previous  models.  In  conjunction  with  the  model  development,  we  arc  collecting 
new  reference  events  with  accurate  locations,  among  them  earthquakes  on  mid-ocean  ridges  and  transforms, 
using  a  constrained  inversion  technique. 
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Objective 

The  aim  of  this  project  is  to  ultimately  improve  the  locations  of  earthquakes  and  other  seismically 
recorded  events  in  order  to  enhance  the  ability  to  monitor  the  Comprehensive  Nuclear-Test-Ban-Treaty 
(CTBT).  Our  strategy  is  based  on  developing  new,  detailed  3-D  models  of  the  mantle,  with  an  emphasis 
on  P  wave  structure.  This  involves  the  construction  of  global  high-resolution  models  with  more  detail  in 
certain  areas  where  particularly  good  data  coverage  is  available.  A  second,  subsidiary  objective  concerns 
the  development  of  additional  techniques  used  in  locating  events  teleseismically  and  regionally  with  sparse 
datasets,  and  with  assessing  the  improvement  in  accuracy  afforded  by  the  new  techniques  and  models.  We 
have  previously  studied  the  location  accuracy  obtained  by  some  previously  existing  mantle  P  wave  models 
[, Antolik  et  al. ,  2000].  Here  we  describe  the  development  of  a  new  joint  P  and  S  wave  velocity  model  of  the 
mantle,  using  a  diverse  and  improved  data  set  of  travel  times  and  surface  wave  phase  velocity  measurements. 

Research  Accomplished 

The  accuracy  of  event  location  with  respect  to  ground  truth  is  not  a  simple  function  of  the  number  of 
free  parameters  or  “resolution”  of  a  particular  velocity  model.  This  was  reported  by  Antolik  et  al.  [2000], 
where  they  tested  four  existing  P  wave  models.  Two  of  these  models,  S&P12/WM13  (hereafter  referred  to 
as  SP12)  [ Su  and  Dziewonski,  1993]  and  MK12WM13  [Su  and  Dziewonski ,  1997]  were  parameterized  in 
terms  of  relatively  low-order  spherical  harmonic  functions.  The  other  two  models,  that  of  van  der  Hilst  et 


Maximum  layer  number 


Figure  1 :  Root  mean  square  mislocation  obtained  from  locating  a  set  of  nuclear  explosions  with  known  locations 
using  model  corrections  from  BDP98  corresponding  to  layers  only  shallower  than  a  certain  depth.  Each  layer  is 
slightly  less  than  200  km  thick  and  are  numbered  in  ascending  order  from  the  top  to  the  bottom  of  the  mantle.  The 
depths  of  670  km  and  1000  km  are  indicated  by  the  vertical  lines.  The  data  point  for  layer  0  is  equivalent  to  locating 
the  events  in  the  reference  1-D  model  PREM  (no  layers  from  BDP98  were  used).  The  majority  of  the  improvement  in 
location  accuracy  is  obtained  by  considering  only  those  layers  above  1000  km. 


al.  [1997]  (HWE97)  and  Boschi  and  Dziewonski  [1999]  (BDP98),  were  parameterized  in  terms  of  constant- 


velocity  blocks  with  horizontal  dimensions  of  2°  and  5°,  respectively,  and  contain  nearly  an  order  of  mag¬ 
nitude  more  free  parameters.  Surprisingly,  the  best  locations  were  provided  by  model  SP12  even  though  the 
models  with  higher  resolution  provide  a  better  fit  to  the  travel  time  residuals.  These  tests  were  performed  us¬ 
ing  a  well-calibrated  set  of  explosions  and  earthquakes  [  Kennett  and  Engdahl,  1991]  and  datasets  consisting 
of  a  limited  number  of  teleseismic  phases  in  order  to  simulate  global  location  of  small  events  by  the  Inter¬ 
national  Monitoring  System  (IMS)  seismic  network.  Model  SP12  provided  consistently  smaller  magnitude 
mislocations  than  the  other  three  models. 

There  are  several  possible  reasons  for  this  result.  The  parameterization  most  often  chosen  for  newer, 
high-resolution  models  is  typically  that  of  blocks  with  constant  velocities  [e.g.,  Boschi  and  Dziewonski, 
1999;  Vasco  and  Johnson,  1998;  Grand  et  al.,  1997].  These  models  therefore  introduce  artificial  lateral 
discontinuities  into  the  structure  and  may  induce  an  unrealistic  shape  in  long-wavelength  anomalies.  A 
lack  of  correlation  between  new,  high-resolution  Earth  models  and  earlier  longer  wavelength  models  has 
previously  been  noted  [Grand  et  al.,  1997],  Another  issue  is  the  coverage  of  the  data  sets  used  in  the 
inversion  problem.  Because  of  the  very  large  number  of  parameters  used  to  define  the  higher  resolution 
models,  most  of  these  models  make  use  of  a  somewhat  limited  data  set  (for  example,  only  phase  travel 
times)  in  order  to  reduce  the  size  of  the  necessary  computer  resources  (memory  and  CPU  time).  This 
may  result  in  a  relative  lack  of  resolution  in  certain  areas  of  the  mantle  (particularly  at  shallow  depths)  to 
which  teleseismic  phase  travel  times  are  not  sensitive.  As  an  example,  we  show  in  Figure  1  the  mislocation 
obtained  from  a  group  of  explosions  using  only  a  subset  of  the  layers  of  model  BDP98  as  a  correction  to 
PREM.  Nearly  all  of  the  improvement  in  the  locations  is  obtained  by  considering  depths  above  1000  km, 
which  shows  that  resolution  of  the  shallow  mantle  is  most  critical  to  teleseismic  event  location.  Models 
which  do  not  recover  the  upper  mantle  structure  well  will  not  provide  the  most  accurate  locations. 
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Figure  2:  Diagram  showing  the  average  amplitide  of  the  diagonal  elements  of  an  inner  product  matrix  for  an  inversion 
for  shear  wave  velocity  [Gu  et  al .,  2000].  Each  cell  corresponds  to  that  portion  of  the  matrix  describing  sensitivity 
of  a  particular  data  set  used  (listed  at  the  top)  to  one  of  the  radial  B- splines  of  the  model.  The  splines  are  numbered 
in  ascending  order  from  the  CMB  to  the  top  of  the  mantle  (i.e.,  spline  1  has  maximum  amplitude  near  the  CMB,  see 
Figure  4).  A  higher  normalized  amplitude  indicates  greater  sensitivity  of  a  dataset  to  the  coefficients  of  a  radial  spline. 


(a) 


(b) 


Figure  3:  (a)  Horizontal  parameterization  of  the  global  model.  The  Earth’s  surface  is  divided  into  362  equal-area 
triangles.  Each  knot  location  corresponds  to  the  center  of  a  spherical  B-spline.  (b)  A  3-D  view  of  a  spherical  B-spline 
used  in  the  lateral  expansion.  The  spline  function  is  localized,  axially  symmetric,  and  has  been  normalized  to  1 . 


It  is  therefore  important  to  use  data  in  tomographic  inversions  which  have  maximum  sensitivity  to  upper 
mantle  structure.  Figure  2  shows  the  average  value  of  the  diagonal  elements  of  an  inner  product  matrix  used 
in  an  inversion  for  S  wave  velocity  [Gu  et  ah,  2000].  Maximum  sensitivity  to  the  upper  mantle  is  provided 
by  the  use  of  surface  wave  dispersion  measurements  and  body  and  surface  waveforms  with  periods  longer 
than  45  s.  Differential  and  direct  travel  times  are  mostly  sensitive  to  lower  mantle  structure.  For  a  P  wave 
inversion  Figure  2  would  be  similar,  without  the  Love  wave  measurements  and  with  the  travel  times  for 
corresponding  P  phases  substituted.  In  deriving  our  new  global  mantle  model,  we  have  therefore  made  use 
of  a  variety  of  datasets. 

For  the  model  parameterization  we  choose  a  combination  of  spherical  and  radial  B-splines.  Horizon¬ 
tally,  the  spherical  splines  arc  centered  at  362  uniformly  distributed  knots  (Figure  3).  The  number  of  free 
parameters  is  close  to  that  involved  in  a  degree- 18  spherical  harmonic  expansion.  Perturbations  to  PREM  in 
both  P  and  S  velocity  arc  represented  by 

=  E  casi  t)  ’  Bi M  (D 

Vq  .  . 

h3 

where  r,  6 ,  and  cj>  are  radius,  colatitude,  and  longitude,  respectively,  Sj  is  the  jth  spherical  spline,  and  Bj  the 
zth  radial  spline.  Cr  j  are  the  model  coefficients  to  be  determined.  The  mathematical  form  of  the  Sj  is  given 
in  Wang  and  Dahlen  [1995].  A  3-D  view  of  one  of  the  spherical  spline  functions  is  shown  in  Figure  3b. 
The  maximum  amplitude  is  localized  to  the  vicinity  of  the  central  knot.  The  spline  approach  thus  combines 


some  of  the  advantages  of  using  a  local,  block  parameterization  with  the  advantage  of  producing  a  smooth 
model. 

Radially,  the  parameterization  consists  of  14  cubic  B -splines  distributed  non-uniformly  with  depth 
(Figure  4).  The  local  support  of  B-splines  allows  for  fast  computation  since  only  neighboring  splines  arc 
needed  to  compute  the  value  of  the  model  at  a  given  point.  The  second  derivatives  of  the  splines  vanish  at 
both  ends.  We  use  the  split  parameterization  of  Gu  et  al.  [2000]  in  which  the  radial  splines  arc  split  at  the 
boundary  between  the  upper  and  lower  mantle  (670  km)  in  order  to  detect  possible  sharp  changes  in  velocity 
variations  across  this  depth.  Six  B-splines  arc  used  for  the  upper  mantle  portion  and  8  for  the  lower  mantle. 
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Figure  4:  Split  radial  B-splines  used  for  parameterization  of  the  global  model.  Two  sets  of  splines  are  used:  6  for  the 
upper  mantle  and  8  for  the  lower  mantle.  This  parameterization  provides  higher  nominal  resolution  near  the  670-krn 
discontinuity. 

We  are  currently  using  a  large  of  number  of  datasets  in  the  inversions  as  well  as  developing  others  for 
future  use.  The  compressional  velocity  model  is  determined  by  travel  times  of  direct  P  and  core  phases 
( PKP ,  PcP)  to  improve  resolution  in  the  lowermost  mantle.  These  datasets  were  adopted  from  those  of 
Engdahl  et  al.  [1998]  but  have  been  improved  through  relocation  of  the  earthquakes  in  a  3-D  model  (SP12). 
The  travel  times  were  then  formed  into  summary  rays  with  sources  at  the  center  of  2°  x  2°  cells.  Events 
deeper  than  50  km  were  grouped  by  depth  into  bins  with  a  thickness  of  100  km.  Events  whose  epicenter 
moved  more  than  1°  from  the  location  published  by  Engdahl  et  al.  [1998]  were  discarded.  The  total  number 
of  summary  rays  is  626,073  for  P,  215,590  for  PKP ,  and  68,473  for  PcP.  We  also  make  use  of  the  Love 
and  Rayleigh  wave  dispersion  measurements  published  by  Ekstrom  et  al.  [1997]  in  the  period  range  35-150 
s,  which  have  their  maximum  sensitivity  in  the  upper  200-300  km  of  the  mantle.  These  were  previously 
used  by  Ekstrom  and  Dziewonski  [1998]  to  derive  a  3-D  model  of  mantle  shear  wave  velocity  expanded  in 


spherical  harmonics  up  to  degree  20  as  part  of  an  earlier  phase  of  this  project.  For  shear  velocity  inversions, 
we  also  use  the  direct  and  differential  travel  time  datasets  described  in  Su  et  al.  [1994].  They  consist  of 
approximately  45,000  measurements  of  S ,  SS,  ScS ,  SS  -  S,  S  -  SKS ,  ScS  -  S ,  and  SKKS  -  SKS  travel 
times.  Future  models  will  also  include  the  extensive  dataset  of  body  and  surface  waveforms  compiled  over 
the  years  at  Harvard. 

We  solve  the  inverse  problem  by  non-linear,  iterative  least-squares  analysis.  The  equations  are  of  the 

form 

<5do  =  A  <5x  (2) 

where  ix  is  the  perturbation  in  the  model  from  the  previous  iteration,  <5do  is  the  corresponding  unexplained 
portion  of  the  data  vector,  and  A  is  the  kernel  matrix.  For  the  final  iteration,  we  require  that 

|A<5x  —  <5do  | 2  +  A2 1  <5x| 2  +  r]2g2  =  min  (3) 

where  A  and  rj  arc  damping  factors  and  g2  is  the  integrated  squared  gradient  of  the  model: 

o  f  Crmoho  0 

9  =  |  y  (Sv/v0)\2drcm.  (4) 
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The  least-squares  model  solution  at  each  iteration  is  computed  by  the  method  of  Cholesky  factorization 
[Trefethen  and  Bau ,  1997], 

An  initial  P  wave  velocity  model  is  shown  in  Figure  5.  With  respect  to  model  BDP98,  the  new  model 
shows  higher  amplitude  velocity  perturbations  in  the  upper  mantle,  which  may  be  the  result  of  the  3-D 
relocation  of  the  sources  prior  to  inversion.  Use  of  the  surface  wave  measurements  produces  lower  velocities 
throughout  the  eastern  and  central  Pacific  than  are  observed  using  travel  times  alone.  At  the  bottom  of  the 
mantle,  the  new  model  has  slightly  lower  amplitude  than  BDP98  and  much  lower  amplitude  than  model 
SP12,  especially  in  the  southern  hemisphere.  This  may  be  partly  due  to  the  relatively  coarse  separation 
of  the  radial  splines  in  the  lower  mantle  in  this  parameterization.  Since,  as  noted  above,  the  upper  mantle 
is  most  critical  for  accurate  earthquake  location,  the  new  model  may  improve  upon  the  location  accuracy 
provided  by  BDP98  and  other  high-resolution  models.  Our  new  models  will  be  tested  using  the  dataset  of 
calibration  events  compiled  at  Harvard  and  other  institutions. 

Constraining  earthquake  locations  on  mid-ocean  ridges  and  transforms 

In  addition  to  the  construction  of  more  detailed  and  accurate  global  velocity  models,  we  have  focused 
our  efforts  on  building  a  larger  database  of  reference  or  “ground-truth”  events  to  be  used  for  model  calibra¬ 
tion.  So  far,  the  geographical  distribution  of  events  with  locations  known  to  an  accuracy  of  5  km  or  better  is 
extremely  limited.  Many  of  these  events  are  nuclear  explosions  and  are  concentrated  in  only  a  few  source 
regions,  making  comprehensive  testing  of  velocity  models  and  calibration  of  new  seismic  stations  very  dif¬ 
ficult.  We  are  experimenting  with  relocation  of  events  on  remote  plate  boundaries  using  a  technique  which 
constrains  the  location  to  a  point  along  the  plate  boundary.  For  example,  by  limiting  the  process  to  only 
large  events  for  which  a  Centroid  Moment  Tensor  (CMT)  solution  exists,  an  assessment  can  be  made  as  to 
whether  an  event  should  be  located  along  a  mid-ocean  ridge  segment  or  an  adjacent  transform  fault  based  on 
its  focal  mechanism.  We  then  perform  a  transformation  of  coordinates  into  the  system  defined  by  the  pole  of 
rotation  between  the  two  plates,  such  that  the  ridge  or  transform  fault  lies  along  a  parallel  or  meridian.  We 
use  accurate  bathymetry  to  determine  the  precise  location  and  trend  of  the  boundary  [Smith  and  Sandwell, 
1997],  In  the  subsequent  location  process,  the  coordinate  (latitude  or  longitude)  which  corresponds  to  the 
plate  boundary  is  held  fixed.  In  this  way,  we  arc  able  to  develop  a  master  set  of  events  which  can  then  be 
used  to  relocate  other  smaller  events  in  the  same  region. 
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Figure  5:  Model  of  P  wave  velocity  at  various  depths  obtained  in  joint  inversion  of  compressional  and  shear  velocity 
using  travel  times  and  surface  wave  phase  velocities.  Scale  shows  fractional  velocity  perturbation  from  PREM. 


We  have  thus  far  compiled  a  catalog  of  these  events  along  most  of  the  transform  faults  in  the  Atlantic, 
Pacific,  and  Indian  oceans.  Figure  6  shows  the  relocation  of  events  along  the  Romanche  Fracture  Zone 
(RFZ)  in  the  central  Atlantic.  The  top  panel  shows  the  ISC  locations  of  events  in  the  CMT  catalog.  After 
relocation  (middle  panel),  the  strike-slip  earthquakes  associated  with  the  RFZ  lie  along  the  trend  defined 
by  the  bathymetry.  The  ISC  locations  arc  systematically  displaced  to  the  south  relative  to  the  RFZ,  and 
the  apparent  mislocations  arc  as  large  as  40  km.  The  results  obtained  thus  far  arc  promising  but  require 
extensive  testing  to  determine  the  accuracy  of  the  relocations. 

Conclusions  and  Recommendations 

Although  the  current  generation  of  global  velocity  models  provides  substantial  improvement  over  1-D 
velocity  models  for  teleseismic  earthquake  location  [Smith  and  Ekstrom,  1996;  Antolik  et  ah,  2000],  there 
is  still  much  room  for  improvement.  The  consistent  accurate  location  of  small  events  to  within  the  area 
specified  for  on-site  inspection  under  the  CTBT  (1000  km2)  will  likely  require  the  use  of  regional  phases  in 
addition  to  teleseismic  travel  times,  and  also  possibly  consideration  of  anisotropy  or  arrival  angle  data.  To 
this  end,  we  are  developing  new  global  models  which  combine  a  local  parameterization  with  the  advantages 
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Figure  6:  Relocation  of  earthquakes  along  the  Romanche  Fracture  Zone  in  the  central  Atlantic.  The  top  panel  shows 
the  focal  mechanisms  of  large  earthquakes  in  the  region  plotted  at  their  ISC  locations.  After  identification  of  the  strike- 
slip  events  to  be  constrained  to  the  RFZ,  these  events  are  relocated  after  coordinate  transformation  (middle  panel).  The 
bottom  panel  shows  the  change  in  location  of  these  same  events,  multiplied  by  a  factor  of  2.5.  Bathymetry  is  taken 
from  Smith  and  Sandwell  [1997]. 
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of  smoothness  over  large  distances.  This  will  enable  straightforward  replacement  of  portions  of  the  global 
model  by  a  more  detailed  model  covering  a  particular  region.  The  smooth  parameterization  will  allow  simple 
computation  of  3-D  raypaths. 

Testing  of  new  global  and  regional  velocity  models,  as  well  as  calibration  of  newer  IMS  stations,  will 
benefit  from  the  compilation  of  reference  events  in  remote  regions.  Our  database  of  relocated  events  on 
mid-ocean  plate  boundaries  will  also  enable  construction  of  more  accurate  travel  time  and  phase  velocity 
datasets. 
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